Two-body momentum correlations in a weakly interacting one-dimensional Bose gas 
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We analyze the two-body momentum correlation function for a uniform weakly interacting one- 
dimensional Bose gas. We show that the strong positive correlation between opposite momenta, 
expected in a Bose-Einstein condensate with a true long-range order, almost vanishes in a phase- 
fluctuating quasicondensate where the long-range order is destroyed. Using the Luttinger liquid 
approach, we derive an analytic expression for the momentum correlation function in the quasicon- 
densate regime, showing (i) the reduction and broadening of the opposite-momentum correlations 
(compared to the singular behavior in a true condensate) and (ii) an emergence of aniicorrelations 
at small momenta. We also numerically investigate the momentum correlations in the crossover 
between the quasicondensate and the ideal Bose-gas regimes using a classical field approach and 
show how the anticorrelations gradually disappear in the ideal-gas limit. 
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I. INTRODUCTION 

Many-body correlation functions contain valuable in- 
formation about the physics of quantum many-body sys- 
tems and therefore their measurement constitutes an im- 
portant probe of the correlated phases of such systems. 
In recent years, ultracold atom experiments have shown 
that atomic correlations can be accessed via many ex- 
perimental techniques, including high-precision absorp- 
tion [ll-Q or fluorescence imaging , detection of atom 
transits through a high- finesse optical cavity @, single- 
atom detection using multichannel plate detectors |8|- 
IllT ] or scanning electron microscopy techniques [l2l |. and 
the measurement of rates of two-bod y (p hotoassociation) 
fl3| or three-body loss processes flil - llol ]. While the loss- 
rate measurements depend only on local correlations, the 
imaging and atom detection techniques typically depend 
on nonlocal correlations which are embedded in the atom 
number fluctuations in small detection volumes (such as 
image pixels) or in the coincidence counts of time- and 
position-resolved atom detection events. 

The development of these techniques have enabled the 
study of a wide range of phenomena in ultracold atomic 
gases, including the Hanbury Brown-Twiss effect @- 
fiol [l7l - [l9j and higher-order coherences [2(| , phase fluc- 
tuations in quasicondensates [2ll. [22j , superfluid to Mott 
insulator transition 0, [H, H, l23|. isothermal compressibil- 
ity and magnetic susceptibility of Bose and Fermi gases 
IS HIS scale invariance of two-dimensional (2D) sys- 
tems [29|], the phase diagram of the ID Bose gas |3fj| . 
l3lj . entanglement and spin squeezing in two-component 
and double- well systems 32h 3a . su b-Poissonian relative 
atom number statistics [35l l37l l38j. and violation of the 
Cauchy-Schwarz inequality with matter waves (39j . 

From a broad statistical mechanics point of view, 
most of these measurements have so far given access 
to either equilibrium position-space density correlations 



or nonequilibrium momentum-space density correlations. 
In this paper, we address the question of equilibrium 
momentum-space density correlations [40j by focusing on 
the two-body correlation function 

G{k,k') = (Sh k Sn k >) = (n k h k ,) - (h k )(n k >), (1) 

for a weakly interacting uniform ID Bose gas. Here, 
5h k = h k — (fife) is the fluctuation in the population n k 
of the state of momentum hk [see Eqs. <j3j , (J4j) , and (|BJ)]. 

To measure Q{k,k') experimentally, one needs to an- 
alyze atomic density fluctuations in a set of momentum 
distributions. Single-shot momentum distributions of a 
ID Bose gas, realizable by confining the atoms to highly 
anisotropic trapping potentials, can be acquired as fol- 
lows. First, by turning off (or strongly reducing) the 
transverse confinement, one ensures that atom-atom in- 
teractions no longer play any role in the system dynamics. 
The longitudinal momentum distribution is unaffected by 
the turning off since, in ID geometry, the turning-off time 
(which is on the order of the period of the transverse con- 
fining potential) is much smaller than the relevant time 
scales of the longitudinal (axial) motion of the atoms. 
Second, the longitudinal momentum distribution can, in 
principle, be measured using an expansion along the long 
axis, after switching off the longitudinal confinement, or 
by using a recently demonstrated technique of Bose-gas 
focusing HlHll (see also [Mill)- 

In the presence of a true long-range order, the Bogoli- 
ubov theory correctly describes the excitations of a Bose 
condensed gas, predicting strong positive correlations in 
Q(k,k') between opposite momenta, k! = —k, for small 
\k\, as shown in Ref. [4(| (see also |47|). However, true 
long-range order is destroyed by long-wavelength fluctu- 
ations in a ID Bose gas J48(; for a large enough system, 
the gas lies in the so-called quasicondensate regime [49| 
where, while the density fluctuations are suppressed as 
in a true Bosc-Einstcin condensate (BEC), the phase still 
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fluctuates along the cloud. In this paper, we show that, 
when the system size becomes much larger than the phase 
correlation length, the positive correlations between the 
opposite momenta vanish. In the thermodynamic limit 
of an infinite quasicondensate, we find an analytic ex- 
pression for Q(k,k') and show that it develops zones of 
anticorrelation on the (k, k') plane. We also analyze the 
crossover from the quasicondensate to the ideal Bose gas 
regime, using a classical field theory, and show how the 
behaviorof Q{k,k') undergoes a continuous transforma- 
tion between the two limiting regimes. 

The paper is organized as follows. In Sec. [TT] we out- 
line the generalities applicable to two-body momentum 
correlations for the uniform ID Bose gas with contact 
interactions. Section IIIII summarizes the known results 
in the regime of a true condensate. In Sec. IIV1 we show 
that the correlations between opposite momenta, which 
exist in the case of a true BEC, disappear in the quasi- 
condensate regime. Here we first use a simple model of a 
quasicondensate (Sec. I IV A[) . followed by the Luttingcr 
liquid approach (Sec. IIVB[) leading to an exact analytic 
result for the two-body momentum correlation function. 
In Section [V] we describe the momentum correlations in 
the crossover from the quasicondensate up to the ideal 
Bose gas limit, using a classical field method. We discuss 
the experimentally relevant aspects in Section IVI1 and 
conclude with a summary in Section IVIII 

II. GENERALITIES 

We consider a uniform gas of bosons interacting via a 
pairwise (5-function potential in a ID box of length L with 
periodic boundary conditions. In the second-quantized 
form, the Hamiltonian density is 

n = ~^^ + ^^^-^t>> ( 2 ) 

where ip(z) and ip ( z ) are the bosonic field operators, m 
is the mass of the particles, g is the interaction constant, 
and p is the chemical potential. In the grand-canonical 
formalism that we are using, the equilibrium density 
p = (ip'ip) is fixed by p and the temperature T, and the 
total number of particles is given by N = pL. Through- 
out this paper, we restrict ourselves to the weakly inter- 
acting regime, which corresponds to the dimensionlcss 
interaction parameter 7 = mg/h 2 p <C 1. 

The momentum distribution (fik) and its correlation 
function G(k, k') are related to the first- and second-order 
correlation functions of the bosonic fields, 

Gx( Zl ,z 2 ) = Gx( Zl - z 2 ) = (ft( Zl )i>(z 2 )}, (3) 

and 

G 2 (z!, z 2 , Z3, Z4) = (^ t (zi)^(z 2 )-0 t (z 3 )'0(z4)), (4) 
via the Fourier transforms 

(n k ) = j Jj L d Zl dz 2 e- lk ^-^G 1 (z 1 ,z 2 ), (5) 



and 

Q(k, k') = jjjj^ dAz e- lk< - zl - Z2 ^e- lk,( - Z3 - z ^ 

x [G 2 (zi, z 2 ,z 3 ,Z4,) - Gx(zi, z 2 )G x (z 3 ,Zi)] , (6) 

where d z = dz\dz 2 dz 3 dz 4 . In Eq. (J3]), the dependence of 
G\(z\, z 2 ) only on the relative coordinate z\ — z 2 follows 
from the translational invariance of the system. 

Several general statements about the momentum cor- 
relation function Q(k,k') can be made, which are valid 
in any regime of the gas. First, the correlation function 
obeys the following sum rule: 

J2Q(k,k') = {N 2 )-(Nr, (7) 

k,k' 

where N = dxip^ (x)ip(x) is the total particle number 
operator. This implies that, within the canonical ensem- 
ble, one has J^k k' S(k, k') = 0. In the grand canoni- 
cal ensemble, the fluctuation-dissipation theorem, which 
connects the particle number variance with the deriva- 
tive of (N) with respect to the chemical potential p 
gives 

V Q(k, k') = k B T^- = k B TL%. (8) 

Second, Q{k, k') possesses several symmetries. In ther- 
mal equilibrium, the position-space correlation functions 
are invariant by the simultaneous refection symmetry of 
all coordinates Zi — > —z%. This symmetry and the bosonic 
commutation relations between the field operators im- 
ply, for periodic boundary conditions, that G(k,k') is 
symmetric around the axis k' = k and around the axis 
k' = -k. 

Finally, for systems that have correlation lengths much 
smaller than the system size L, the two-body momentum 
correlation function G(k, k 1 ) can be split into a 'singular' 
part and a regular function. (We use the term 'singu- 
lar' in the sense of the Kronccker (5-function, which turns 
into the Driac (5-function singularity in the thermody- 
namic limit of L — > 00.) To show this, let us first note 
that if we assume the existence of a finite correlation 
length 1$ for the decay of the first-order correlation func- 
tion G\ (21,2:2), then the second-order correlation func- 
tion G 2 (zi, 22, 23, 24) must have the following two asymp- 
totic limits: 

G 2 (2l, 2 2 , 2 3 , 24) ~ Gl(2i — 2 2 )Gl(2 3 — 2 4 ), 

(9) 

for \zi - z 3 \ > 1$ and \zi - z 2 \, \z 3 - z 4 \ < 1$, 
and 

G 2 (z l7 Z 2 ,Z 3 ,Z 4 ) ~ Gl(2l ™ 2 4 )Gl(2 2 - z 3 ) 

+G 1 (z 1 -z A )S(z 2 -z 3 ), (10) 
^ for \z\ - z 2 \ > 1$ and \z\ - 2 4 |, \z 2 - z 3 \ < 1$. 
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In Eq. (fTO|) . the 6- function term appears simply as a 
result of normal ordering of the operators in Eq. (QJ. 
By separating out the two asymptotic limits, Eqs. ([9]) 
and (jlOp , we can write 

G 2 (zi, z 2 , z 3 , Zi) = Gi(zi - z 2 )Gi(z 3 - Z4) 
+ G 1 (z 1 - z 4 )G 1 (z 2 - z 3 ) 
+ Gi(z 1 - z 4 )6(z 2 - z 3 ) 
+ G 2 (zi, z 2 ,z 3 ,za), (11) 

where G 2 (z\, z 2 , Z3,z 4 ) is the remainder term. 
By substituting Eq. (fl~Tj) into Eq. ^ we obtain 

Q(k, k') = ((h k ) + (n k ) 2 )6 k . k , + G(k, k'), (12) 

which shows explicitly that Q(k, k') can be written down 
as a sum of a singular and regular contributions. The 
first term in Eq. (fT2j) is the shot noise, the second term 
is the bosonic "bunching" term, which describes the ex- 
change interaction due to Bose quantum statistics, and 
the last, regular term Q(k,k') [the Fourier transform of 
G 2 (z\, z 2l Z3, Z4)] describes the exchange of momenta be- 
tween the particles during the binary elastic scattering 
processes and is nonzero only for an interacting gas. For 
nonintcracting bosons, Wick's theorem can be applied di- 
rectly to the G 2 (z±, z 2 , Z3, Z4) function, which then leads 
to a vanishing G 2 (zi, z 2l z 3l z 4 ) and hence only to the sin- 
gular terms in Eq. (IT21) . 



III. TRUE CONDENSATE 

At T = 0, the first-order correlation function de- 
cays algebraically as G\{zx,z 2 ) ~ (£/|zi — z 2 \)^/ 2 ^ for 
\zi — z 2 \ ^ i [50 - 4521 ] . where £ = h/^/rngp is the healing 
length. As 7 <C 1 in the weakly interacting regime, the 
algebraic decay is very slow, leading to an exponentially 
large phase correlation length. Indeed, defining i^ ' as 
the length for which G\ (zi,z 2 ) decreases by a factor of 
e, we find g| 



(0) 



(13) 



At finite temperatures, the algebraic decay of 
Gi(zi,z 2 ) remains valid for distances £ <C \z\ <C It, 
where It = h 2 ImknTt, = (H/kfiT)\/qp/m is the phonon 
thermal wavelength [5l|, [52j]. For distances \z\ ^> It, 
on the other hand, the correlation function decays ex- 
ponentially [see Eq. (f2~5j) below] with the characteristic 
temperature-dependent phase coherence length 



1<I,{T) = h 2 p/mk B T. 



(14) 



Considering now a system of size L <§; minl^ -*, l^}, we 
can assume true long-range order in the system and use 
the Bogoliubov theory to describe the momentum corre- 
lations as was done in Refs. [3,|43|. We briefly recall the 
relevant results here. 



The momentum correlation function G(k,k') is differ- 
ent from zero only for k = k' and k = — k' . For equal 
momenta k = k' , one finds Q{k,k) — (n k ) + (h k ) 2 , 
which is similar to the ideal Bose gas behaviour, ex- 
cept that the standard Bose occupation numbers (n k ) = 

( e (E k -tl)/k B T _ -J^-l re pl aCeC [ Q y 



(n k ) = (l + 2n fc ) 



E k +gp 
2e fc 



(15) 



Here e k = E k (E k + 2gp) is the energy of the Bogoli- 
ubov modes, E k = h 2 k 2 /2m is the free particle disper- 
sion, and n k = (e £fc / fcsT — l) are the mean occupation 
numbers of Bogoliubov modes. For opposite momenta, 
k = —k', one has 



Q(k,-k) = (l + 2n fc ) s 



9P 
2e fc 



(16) 



A convenient way to characterize the relative strength 
of the opposite and equal momentum correlations is via 
the normalized pair correlation function 



V{k) 



G{k,-k) 
G(k,k) 



1 



({h k - n_ fc ) 2 ) 
2(Sh 2 ) 



(17) 



Here, V(k) = 1 corresponds to perfect (maximum) corre- 
lation between the opposite momenta, whereas V(k) = 
corresponds to the absence of any correlation. 

At T = 0, one has n k = and G(k 1 —k) = G(k, fc), and 
therefore the Bogoliubov theory predicts perfect corre- 
lation between the opposite momenta, V{k) = 1. Such 
perfect correlation stems from the fact that the depletion 
of the condensate in the Bogoliubov vacuum simply cor- 
responds to the creation of pairs of particles with equal 
but opposite momenta. 

At finite temperatures, n k is different from zero, nev- 
ertheless the normalized pair correlation is still close to 
its maximum (perfect correlation) value, V(k) ~ 1, for 
phonon excitations with k <C l/£ for any value of n k . 
This can be understood from the fact that the phonons 
are mainly phase fluctuations, so that they correspond to 
equal- weighted sidebands at momenta k and — k of the 
excitation spectrum. On the other hand, for particle- 
like excitations, with k 3> l/£, the thermal population 
of particles leads to a decrease of 'P(k). More precisely, 
for l/£ < k < y/rnksT /h, which corresponds to particle- 
like excitations whose occupation numbers are large, one 
obtains 'P(fc) <C 1. Finally, at very large momenta, 
k 3> ^/mk~BT/h, for which the occupation numbers are 
negligibly small, one again recovers the zero-temperature 
result P(k) ~ 1. 



IV. QUASICONDENSATE REGIME 

A. Effect of phase fluctuations 

The above results obtained using the Bogoliubov the- 
ory are valid when the temperature is small enough so 
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that the phase correlation length is much larger than the 
system size, 1^ 3> L. While this condition is easier to sat- 
isfy in 3D or quasi- ID systems, it is generally not fulfilled 
for purely ID gases. 

In a large enough ID system or at high enough tem- 
peratures, the long-range order is destroyed by long- 
wavelength phase fluctuations, having a characteristic 
temperature-dependent correlation length l^. When 
l<f, -^l L, such a system is said to enter into the so-called 
quasicondensate regime [49j , in which the density fluctu- 
ations are suppressed while the phase still fluctuates. As 
we show here, the two-body correlation between opposite 
momenta is expected to vanish in the quasicondensate 
regime. 

To give a crude, yet simple estimate of the two-body 
momentum correlations, we can divide the system into 
domains of length 1$ and assume that (i) within each do- 
main, the spatial variation of the phase is small, and 
therefore the Bogoliubov approach for a true conden- 
sate can be applied to each domain, and (ii) the relative 
phases between two different domains are uncorrelated. 
For each domain, indexed by a, the field operator ip a (z) 
can be expanded according to the Bogoliubov theory, 

M*) = e # ° f VP + 4= E 5 ^,k e~ lkz I , (18) 

where the first term is the mean- field component, the 
second term is the fluctuating component expanded in 
terms of plane- wave momentum modes Sip a ^, <Pa is the 
mean global phase of the domain assumed to be a random 
variable distributed uniformly between and 27r, and the 
summation is over the momenta that are quantized in 
units of 27r//0. 

Using the fact that the momentum component tpk = 
Jq dzip(z)e lkz of the full field -tp(z) can be decom- 
posed as ipk = \/l<f>/ LJ2 a ^a,k& %l ^ a for k ^ 0, we obtain 
the following expression for the momentum correlation 
function: 

(n fe rW) = ( l A y, Wl, A*'> 

^ ' a/3jS 

x (k,k'^0). (19) 

Here, the overline above the exponential factor stands 
for averaging over the random mean phases of different 
domains. 

Within the Bogoliubov theory, the Hamiltonian is 
quadratic in 5ip at k and one can use Wick's theorem 
to evaluate the four-operator correlation function in 
Eq. (|T9"|) . Only pairs of operators belonging to the same 
domain give a nonzero contribution since different do- 
mains are uncorrelated. Among these pairs, only terms 

( Sq Pi,k^a,k)' (fya,k5i>a,-k), and {S^ k Si>l,_ k ) survive. 

To evaluate these terms in the most transparent way 
we make use of the classical field approximation [53| (see 



also Sec. IV Ap . treating the operators 5tp ay k and 6ffi a k as 
c-numbers, dtpa^k and 8ij)* a fc , and assuming that the re- 
spective mode occupations (hk) = (<5^* t.6i/) a k ) are mu ch 
larger than one, (Hk) ~> 1. For fe<Cl/£, the excitations in 
each domain are almost purely phase fluctuations so that 
6ip a - k = ~Sip* ak and therefore (Sip^kWa-k) - ~( h k) 
As a result, for the regular part of the momentum 
correlation function we obtain 

G(k, k') ~ 6 k ,-k- 2 £ e-^~M(n k ) 2 . (20) 

Averaging over the phases gives e _2l (* Q_< ^) = d a ^, 
which singles out only the diagonal in a and j3 terms; 
there are L/l$ such terms in the sum in Eq. (|20[) . Ac- 
cordingly, for the correlation function with opposite mo- 
menta, we find Q(k, —k) ~ (l^,/L)(nk) 2 , whereas the cor- 
relation function for equal momenta is given by Q(k 1 k) = 
(nk) 2 + (rik) — (nk) 2 , for (n k ) > 1. Therefore, for the 
normalized pair correlation V(k) we obtain the following 
simple result 

V(k) ~ /£ ^<1, (21) 

which shows that the correlations between the opposite 
momenta are inversely proportional to the system size L 
and therefore are vanishingly small for L 1$. 

The above simple model is not capable of capturing 
features of G(k,k / ) on momentum scales smaller than, 
or of the order of, the inverse phase correlation length, 
k < l/l<p. For such momenta, the two-body correlation 
function is calculated below using a more rigorous Lut- 
tinger liquid approach. The results obtained within this 
approach confirm the simple scaling behaviorobtained in 
Eq. (|2Tj) . Moreover, the Luttingcr liquid results show 
that the correlation function between different momenta 
is no longer singular on the antidiagonal k! = ~k and 
that it develops zones of anticorrclation. 

B. Two-body correlations in the Luttinger liquid 
approach 

The condition for the quasicondensate regime [5|| is 

In this regime, the correlation functions in Eqs. ([3]) 
and (j?]) arc dominated by the long-wavelength (low- 
energy) excitations and the Hamiltonian density reduces 
to that of the Luttinger liquid [lH, [56| : 

n L = 9 -(5p) 2 + ^(d z $) 2 . (23) 

Here, 5p(z) is the operator describing the density fluctu- 
ations, canonically conjugate to the phase operator <j>(z), 
with the commutator [Sp(z), 4>(z')\ = iS(z — z'). 
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The density fluctuations are small in the quasiconden- 
satc regime and, as long as the relative distances con- 
sidered are much larger than the healing length £, they 
can be neglected when calculating the correlation func- 
tions © and (gl) [H[5l[53,[5i|. As the Luttinger liquid 
Hamiltonian given by Eq. (f2"B"|) is quadratic in 0, the first- 
order correlation function G x {z l ,z 2 ) = p{e i ^ z ^~^ z ^) 
can be expressed through the mean-square fluctuations 
of the phase using Wick's theorem: 



G 1 (z 1 ,z 2 )=pe-i«*^-^» 2 >. 



(24) 



Neglecting the contribution of vacuum fluctuations 
compared to thermal ones [H^ the calculation of the 
mean-square phase fluctuations leads to an exponentially 
decaying first-order correlation function (5lL l52l [58| , 

G 1 (z 1 ,z 2 )=pe-I zi - Z2 l/ 2 ^, (|z 1 -z 2 |»0. (25) 

This defines the finite-temperature phase coherence 
length 1$, given by Eq. (|T4l . and leads to a Lorentzian 
distribution for the momentum mode occupation num- 
bers, 



1 + {2l^kf 



(26) 



valid for k -C l/£. 

Similarly, the two-body correlation function 

G 2 (z 1 ,z 2 ,z 3 ,z 4 ) = p>( e i[$(zi)-foz2)+4>(z3)-$M]) can 
be represented in terms of the first-order correlation as 



G 2 (Z!,Z 2 ,Z 3 ,Z 4 ) 

_ Gi(z!- z 2 )Gi(z 3 - z 4 )Gi(zi- Zj)Gi(z 2 - z 3 ) 
Gi{zi - z 3 )Gi(z 2 - Zi) 



(27) 



By substituting Eqs. (J25j) and j2j]) into Eq. ©, we 
find that the two-body momentum correlation function 
indeed has the form of Eq. (fT2|) [60l ] , in which the regular 
part G(k, k') can be written as 

§(k, k') = l j- {pl^f J-(2^fc, 21^') , (28) 

where F(q, q') is a dimensionless function given by 

T{ i \ = 256 

[9,q ' (q2 + iy(q>2 + m(q + q>)2 + 16} 

X [(q 2 + 5qq' + q' 2 )qq' - 2(q 2 - qq' + q' 2 ) - 7] , (29) 

with q = 2/^fc and q' = 2/^fc'. We note that the re- 
striction of these results to k -^i l/£ implies q <C — 
2T co /T, and that the scaling of G(k, k') with the inverse 
size of the system L coincides with the one obtained in 
Eq. (HU). 

We now wish to check the constraints on the function 
Q{k,k') imposed by the sum rule, given by Eq. ((5J. In 
evaluating the different terms in the left and right hand 
sides of Eq. © , we note that (i) for the derivative term 
we can use the equation of state for the quasicondensate 



regime, p = fi/g, (ii) the term ^2 k (fik) [coming from the 
singular part of Q(k,k')] is given simply by ^fc(^fc) = 
N = pL, and (iii) the term *^2 k {h k ) 2 can be evaluated 
using Plancherel's theorem and Eq. (|2"5"]) . As a result, 
the sum rule is reduced to I6ll 



(27T) 



^2 Jj dqdq'J : {q,q') 



T 

± c.c. 



(30) 



where the contribution of the X}fc(^*) * s ignored on the 
grounds that it is of the order of T/T^ (where = 
H 2 p 2 /2mfcs), which is always much smaller than unity 
in the entire range of temperatures T < T co . 

Th evaluation of the integral on the left-hand side of 
Eq. (|30|) gives the value of —8, implying that the sum 
rule is indeed approximately satisfied as long as T <C T co , 
i.e., deep in the quasicondensate regime. On the other 
hand, as the temperature increases and approaches the 
quasicondensation crossover T co , the term (T/T co ) 2 on 
the right-hand side of Eq. (|3"U)) becomes non-negligible, 
implying that our result for the pair-correlation function 
Q(k,k'), given by Eqs. ([28)) and (I29|) . is no longer valid 
as it fails to satisfy the sum rule [62( . The physical origin 
of this failure lies in the fact that the density fluctuations 
at temperatures near T co are no longer negligible. 

The two-body correlation function Q(k,k'), given by 
Eq. (TT2")) . in the quasicondensate regime, of which the 
regular part G(k 1 k / ) is described by the universal di- 
mensionless function F(q, q'), is one of the key results of 
this paper. The function !F{q 1 q') is shown in Figs. [Ha)- 
(c); as J-(q,q') is independent of the system size L, it 
essentially describes the (unnormalized) two-body mo- 
mentum correlations in the thermodynamic limit. As we 
see, the correlation function is nonzero on the entire 2D 
plane of momentum pairs (k, k'); this can be contrasted 
with the singular behaviorof correlations in the true con- 
densate where Q{k,k') was nonzero only for k' = ±k. 
This effective broadening of correlations is the first conse- 
quence of large phase fluctuations in the quasicondensate 
regime compared to the behaviorin the true condensate. 
Next, Q(k,k') and V(k) both scale as l^/L and there- 
fore are vanishingly small as I > For k = —k', this 
means that the perfect opposite-momentum correlations 
\P{k) = 1], which were present in the true condensate, 
essentially disappear in the quasicondensate regime. Fi- 
nally, we find negative correlations (or anticorrelations) 
in Q(k, k'); these are pronounced mostly in the regions of 
k'k < [see Fig. QJb)]. The correlations fall to zero on 
a typical scale of k ~ 1/1$ (q = 2l^k ~ 1). This is ex- 
pected, as l<f, is the length scale governing the first-order 
spatial correlation function G\{zx,z 2 ) in the quasicon- 
densate regime, and the momentum correlations depend 
only on G±(zi, z 2 ) in this regime. 

To gain further insights into the strength of the two- 
body correlations, we consider the normalized regular 
part of the two-body correlation function, 



s (2) (M') 



G{k,k') 

{h k )(h k ,) 



(31) 
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Figure 1. (Color online) (a) Dimensionless regular part of 
the (unnormalized) two-body momentum correlation func- 
tion, T(q,q'), given by Eq. (|29|) . of a uniform ID Bose gas in 
the quasicondensate regime, (b) Same as in (a), but showing 
the details at small correlation amplitudes (see the scale on 
the color bar) and in a larger window of values of (q, q'). The 
small negative and positive amplitudes seen here get "mag- 
nified" when the function F^q,^) is normalized to {fik){n' k ), 
as is done in Fig. [2] (c) Function J-(q, q') along the diagonal 
(q = q 1 ) and antidiagonal (q = —q'). 



Using Eqs. I|28p and (J2S|), this can be rewritten as 



(32) 



V 




Figure 2. (Color online) Normalized regular part of the two- 
body correlation function f(q,q') as a function of the dimen- 
sionless momenta q — ll^k and q' = H^k! . 



where 



•F(g,g') 
16 



1 + I 2 ) 1 



(33) 



is a dimensionless universal function describing the two- 
body correlations of a ID quasicondensate in the thermo- 
dynamic limit. The function f(q,q') is plotted in Fig [5] 
As we see, the normalization leads, at k 3> 1/Z^,, to the 
recovery [cf. Fig. QJb)] of positive correlations around the 
the antidiagonal k' = —k, predicted by the simple model 
of Sec. IIV Al These correlations can be also thought of as 
the remnants of the nearly perfect correlations in a true 
condensate at k -C 

Finally, we note that g^ 2 '(k,k') can be related 
to Glauber's normally-ordered second-order correlation 
function 



2 (2) (M') 



(34) 



Indeed, by reordering the creation and annihilation oper- 
ators, we can first express the g^ (k, k') function in terms 
of the correlation function G(k, k'), given by Eq. ([T]): 



g^(k,k') = l- 



Sk,k' + 



G(k,k') 
(h k )(n k >) 



(35) 



Using now the general structure of Q{k, k') from Eq. (fT2 
(valid for L 3> 1$) we obtain 



g^\k,k') = l + 5 k . kl +g^(k,k') 



(2)/ 



(36) 



Here, the first term corresponds to uncorrelated atoms, 
the second term is the bosonic bunching term, and the 
last term is the normalized regular part corresponding to 
g{k,k') given by Eq. 

According to our results, the normally ordered nor- 
malized correlation function for equal momenta is given 
by 



g^(k, k) = 2 + gW(k, k) = 2 + 0(l+/L), 



-,(2) 



(37) 
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while for opposite momenta it is given by 

g {2) (k,-k) = l + g {2} (k,-k) = l + Ofo/L). (38) 

The small contributions 0(L,/£) are described by Eq. 
(|32[) and are, in principle, detectable using the precision 
of currently available experimental techniques. Apart 
from the need for high precision on the signal, resolv- 
ing the shape of the g( 2 '(k, fc')-function requires experi- 
mental momentum resolution better than the separation 
Afc = 2tt / L between the individual momentum states 
(for resolutions that are insufficient to resolve separations 
of - 1/L, see Sec. ED- 

As we see from Eq. (|3T|) . the amplitude of equal- 
momentum correlations is close to the pure ther- 
mal bunching level of g( 2 >(k,k) = 2, implying large 
momentum-space density fluctuations. The nearly ther- 
mal level of correlations here is due to the large phase 
fluctuations present in a ID quasicondensate. This makes 
the equal-momentum correlations analogous to those of 
a 'speckle' pattern fl8j where many sources with random 
phases contribute to the familiar Hanbury Brown-Twiss 
interference We emphasize, however, that the nearly 
thermal equal-momentum correlations are obtained here 
for a quasicondensate, which should be contrasted to 
the uncorrelated level of the two-point correlation func- 
tion in position space [5^, |57|, g^(z,z) — li due to 
the suppressed real-space density fluctuations. Equa- 
tion (f3"5)) . on the other hand, shows that the opposite- 
momentum correlations arc close to the uncorrelated 
level of g^ (k,—k) ~ 1, which is in contrast to the strong 
respective correlations [g^(k, —k) = 2+l/(n/-) at T = 0, 
and g^(k, —k) ~ 2 at finite T for (hk) 3> 1] present in a 
true condensate. As we mentioned earlier, the opposite- 
momentum correlations arc essentially destroyed by the 
phase fluctuations. Finally, a significant region of pairs 
of momenta k' ^ k around the origin shows a small de- 
gree of anticorrelation, g^ 2 >(k,k') < 1, which was not a 
priori expected. 



V. FROM THE QUASICONDENSATE TO THE 
IDEAL BOSE GAS REGIME 



A. Classical field approach 

In the quasicondensate regime, T <C T co , higher-order 
correlation functions can always be expressed in terms of 
the first-order correlation function as in Eq. ([2~T|) . There- 
fore, G(k, k') in Eq. ([5]) contains the same information 
as the momentum distribution (rife), given by Eq. (JSJ. 
In particular, the dependence on the temperature comes 
about only through the phase correlation length 1$. This 
is, however, no longer true when the temperature be- 
comes of the order of the crossover temperature T co , in 
which case the physics depends not only on the phase 
fluctuations, but also on the density fluctuations. 

To compute the correlation functions at T > T co , we 
resort to the classical field (or c-ficld) approach of Rcf. 



[53| . In this approach, the quantum field operators ip and 
^ are approximated by c- number fields ip and ip* i whose 
grand-canonical partition function (in a path-integral for- 
mulation) is given by 



Z = \Vi\)Vif exp - 



1 



k B T 



dz H c 



(39) 



Here the function H c (ip, -0*) is obtained from the Hamil- 
tonian density ([2]) by replacing the operators with c- 
fields. The classical field approach is expected to be valid 
for high occupancy of the low-momentum modes con- 
tributing to the momentum correlation function. This 
condition is satisfied in a broad range of temperatures, 
including in the quasicondensate regime, gpe~ 27! / V 5 ? < 
ksT < ^/ 7 )1i 2 p 2 /ni [Hsl, [63|, and up to the temperatures 
corresponding to the degenerate ideal Bose gas regime, 
y /jh 2 p 2 /m < k B T < h 2 p 2 /m O- 

It is convenient to introduce a dimensionless field ip = 
i^/tpo an d a dimensionless coordinate s = z/zq , with 



4>o 



mk 2 B T 2 



1/6 



m 2 gkBT 



1/3 



(40) 



and rewrite the effective 'action' in Eq. (|39| in the di- 
mensionless form: 



— j Q dzU c = J ds ^-\dM 2 + 2^! 4 " ^l 2 



(41) 



This form of the action is controlled by a single dimen- 
sionless parameter 



n 



mg 2 k 2 B T 2 



1/3 



/'■ 



(42) 



Because of the scaling relations (|40|) , the density p = 
{^l>*ip) can be written as p = h(j])(mk 2 B T 2 jY?g) x ^ using 
a dimensionless function h(rj) = ( - 0* - 0). Similarly, the 
phase correlation length 1$, given by Eq. (jT4j) , can be 
written as 1$ = z h(r]). Thus, the length scale Zq can 
be replaced by and therefore the one- and two-body 
correlation functions in Eqs. © and (@| scale as 



G 1 (z l ,z 2 ) = phx ( -i, -^517 
tit, tit, 



(43) 



n 1 \ 2 > i z i z 2 z 3 z 4 \ 

&2[Zi,Z2,z 3 ,Z4) = p h 2 (— ;ti) , (44) 

\i(j> l(j> l(j> i(j> J 

where h\ and h 2 arc dimensionless functions. 

By substituting these scaled correlation functions into 
Eq. © and using the fact that the integrand is invariant 
by a global translation of the coordinates, we find 

Q(k, k') = l j-(pl 4> ) 2 T(2l^ 21^'; rf), (45) 



8 



where J-(q, g'; rf) is a dimensionlcss function parametrized 
by 77. This relation generalizes Eq. ([25)1 beyond the qua- 
sicondensate regime, with the departure being character- 
ized by the value of 77 (see below). 

To find T(q, q'\ 77), we still need to calculate the di- 
mensionlcss function h and correlations h\ and /12 for 
the rescaled fields ijj and -0* with the action given by 
Eq. (|41j). As shown in Ref. [53|, this c-field problem can 
be mapped into the quantum-mechanical problem of a 
particle moving in an external potential. More precisely, 
expressing the action (|41j) in terms of the real and imag- 
inary components of ip = x + iy and interpreting s as 
the imaginary time, the problem can be mapped to the 
quantum mechanics of a particle in two dimensions with 
the Hamiltonian 

H=\ & + pI) + \{* 2 + y 2 f- V (x 2 + y 2 ) . (46) 

Calculating the eigenvalues and matrix elements of this 
Hamiltonian allows one to compute the correlation func- 
tion J^(q, q'; 77). This is done in Appendix [Al 



B. Correlations in the crossover region 

The power of the c-field approach lies in the ability 
to describe the momentum correlations not only in the 
quasicondensate regime, but also in the entire crossover 
region between the quasicondensate and the degenerate 
ideal Bose gas. As shown in Appcndix[Bl in the quasicon- 
densate regime where 77 3> 1 (corresponding to a positive 
chemical potential /j), we recover the results of Sec. UVBl 
with Eq. O referring to F{q, q'; +00) = F{q,q'). The 
opposite limit 77 <C — 1 corresponds to the degenerate 
ideal Bose-gas regime with negative /z. In this case, the 
quartic term in the Hamiltonian (|4"6]l has a negligible ef- 
fect on the lowest energy eigenstates and the problem is 
reduced to a simple two-dimensional harmonic oscillator. 
As shown in Appendix [C] in this limit, we obtain the 
ideal Bose-gas result of T(q, q 1 ; —00) = 0. 

In Fig. E)a) we show how the antidiagonal correla- 
tion function J 7 ^, — q; 77) changes from its quasiconden- 
sate value of Eq. ([29)1 to zero as 77 is continuously changed 
from +00 to —00. To quantify the width of the crossover 
in terms of 77, we consider the peak value of the corre- 
lation function, .F(0,0; 77), and plot it as a function of 77 
in Fig. [3lb). As we sec, J-"(0,0;7/) goes from its mini- 
mum value of about —112 in the quasicondensate regime 
(77 1) to zero in the ideal Bose gas regime (77 <C —1). 

We can define the crossover region to correspond to 
771 < 7/ < 772, where the bounds 771 and 772 arc chosen, 
respectively, at 20% and 80% of the value of ^"(0,0; 77) 
in the quasicondensate regime; our numerical solutions 
give then 771 ~ —1.1 and 772 — 2.0. Recalling that 
the dimensionless parameter 77 is defined via Eq. (|42p. 
this can be converted into the crossover bounds on the 
chemical potential, 771 < pjp ca < 772, where p co = 
k B T(mg 2 /h 2 k B T) 1/z [11]. Similarly, recalling that the 




Figure 3. (Color online) (a) The antidiagonal correlation 
function F(q,—q;r)), for (from top to bottom, dashed lines) 
77 = —1.87, 0, 1.12, and 2.56. The lowest (solid) curve is for the 
limiting quasicondensate regime, !F{q,—q) = J-(q,—q;+oo), 
described by Eq. f)2Qp and shown in Fig. [T](c). (b) The mini- 
mum value of F(q, —q; 77) as a function of 77. The shaded area 
shows the crossover region between 771 < 77 < 772 (see text). 



density p was determined by the dimensionless function 
h(rf), via p = h^^mkgT^/hPg) 1 / 3 , we can use the nu- 
merically found values of h(rj) to rewrite the crossover 
bounds in terms of the density as 0.5 < p/p co < 1-6, 



where p co = (mk B T 2 /tfg) 1 / 3 



VI. EXPERIMENTAL CONSIDERATIONS 

The results obtained so far are directly applicable to 
experimentally measured momentum distributions and 
correlation functions as long the momentum resolution 
is sufficient to resolve the individual momentum states 
separated by Ak = 2ir/L. However, typical resolution in 
ultracold-atom experiments is insufficient to resolve mo- 
mentum scales of the order of 1/L. Because of this, the 
measured signal corresponds to integrated atom-number 
counts Nk in individual detection "bins" (such as cam- 
era pixels in absorption imaging) corresponding to the 
momentum k. 

To address the situation with low momentum resolu- 
tion and relate our calculated correlation functions to 
the experimentally accessible quantities, we assume that 
the detection bin size in momentum space fulfills 
Afe ^> 1/L. In addition, we assume that <C 1/^0, 
so that the bulk of the momentum distribution is still 
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well resolved. With these assumptions, the average (over 
many experimental runs) atom number in a bin (N^) is 
related to the original average mode occupation number 
(hk) via (N k ) = (LA/ c /27r)(nfc), i.e., it accounts for a 
factor equal to the number of original momentum states 
contributing to the bin, Afc/Afc = L/S. k /2-K. Next, the 
average correlation between the bin population fluctua- 
tions is related to the correlation function Q(k, fc'), given 
by Eqs. (12]) and ([28]), via 



(N k N k ,) - (N k )(N k ,) = (N k )S k , k 
+ (N k )(N k 



Ll\ k L 



(47) 



where the universal function f(q, q') is given by Eq. (|33l) . 
In Eq. (|4"T|) . the first term is the shot noise, the second 
term corresponds to the bunching term in Eq. (|12[) , and 
the third term is the contribution of the regular part, 
G{k,k'). 

The shot-noise term is much smaller than the bunching 
term as long as highly populated momentum states are 
considered, i.e., (h k ) S> 1. The latter condition is satis- 
fied for momenta k < 1//^ (containing the bulk of the mo- 
mentum distribution), and therefore the shot-noise term 
can be safely neglected for these momenta. 

Comparing now the bunching term and the regular 
component (with the comparison being relevant only for 
k = fc'), we see that they both scale inversely proportion- 
ally to the system size L (and therefore are small), but 
the regular component is much smaller than the bunch- 
ing term as the dimensionlcss function /(2/^fc, 21^') is 
of the order of one and we have assumed Afc <C l/l^. 
However, the ratio of these two terms is independent of 
L and therefore is finite in the thermodynamic limit. As 
this ratio is proportional to A^Z^ <C 1, detecting the 
contribution of the regular component is going to de- 
pend on actual experimental parameters and the preci- 
sion (signal-to-noise) with which the atom-number fluc- 
tuations can be measured. High-precision measurements 
of atom-number fluctuations, capable of resolving small 
signals like this or even below the shot-noise level, have 
been demonstrated in many ultracold atom experiments 

m him m m m. 

Considering now the cross correlation between atom- 
number counts in different bins, k' ^ k, we see that 
the only contribution to Eq. (j4"T|) comes from the reg- 
ular component. This scales as 1/L, but again such a 
magnitude of the cross correlation should be accessible 
with state-of-the-art measurement techniques, as demon- 
strated, e.g., in Ref. [24| . 

VII. SUMMARY 

To summarize, we have calculated the two-body mo- 
mentum correlations for a weakly interacting, uniform 
ID Bose gas. Our results span the entire quasicon- 
densate regime, where the correlations are derived an- 
alytically in terms of a universal dimensionlcss function 



/(2^fc, 21^'), as well as the crossover to the ideal Bosc- 
gas regime, where the correlations are calculated numeri- 
cally using the classical field method. A natural extension 
of the approaches employed here would be the calculation 
of these correlations for harmonically trapped gases [66[ , 
which is more appropriate for quantitative comparisons 
with experiments beyond the widely used local-density 
approximation. Calculating and understanding the mo- 
mentum correlations in the strongly interacting regimes 
would require the development of alternative theoretical 
approaches and remains an open problem. The knowl- 
edge of such correlations is important in the studies of 
nonequilibrium dynamics from a known initial state and 
the subsequent thermalization in isolated quantum sys- 
tems [U(23. 
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Appendix A: Classical field approach: diagonalizing 
the effective hamiltonian 



In this appendix, we outline how the classical field cor- 
relation functions can be computed using the the equiva- 
lent quantum-mechanical problem of a particle in an ex- 
ternal potential. We recall that the classical-to-quantum 
mapping is done by expressing the c-field ip = x + iy via 
its real and imaginary parts which, in turn, are treated 
as coordinates of a quantum-mechanical particle in imag- 
inary time with the Hamiltonian given by Eq. (|46p . Here 
we imply the scaling of Eq. (|4T)]) but omit the tilde on top 
of the coordinates and c-fields for the sake of notational 
simplicity. 

Using the notations of effective quantum mechanics, 
the first- and second-order correlation functions of the 
c-ficlds, 



(Al) 



and 



G 2 ( Sl , S 2 , S 3 , 8t) = (^* ( Sl ) (s 2 ) r (s 3 ) V (84.)) , (A2) 

are given by 



Gi 



Tr[U L 



(A3) 
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and 



Go 



TrlU L _ s ,y 1 U s , i _ s ,y 2 U s ,_ s ,y 3 U 4 _ s ,* 4 U s , i ] 
Tr[U L ] 



(A4) 

where we have omitted the arguments of G\ and G 2 for 
notational brevity. Here U s = e~ sH is the imaginary- 
time evolution operator generated by the Hamiltonian 
Pr?|) . In Eq. (|A4[) . we take into account the automatic 
time ordering implied by the path integral by introduc- 
ing s[ > ... > s 4 , which is the ordered permutation of 
81, ...,84. The operator stands for ip = x + iy if s' fc 
equals s 2 or S4, or for ip* = x — iy if s' k equals si or S3. 

In the limit L — > 00, the ground state |0) gives the 
dominant contribution to both the numerator and de- 
nominator in Eq. (|A4[) , and hence the correlation func- 
tions reduce to 



Gi = (o|*itr ai _ 8i * 2 |o) 



(A5) 



and 



G 2 = (0|*if/ si _ 4 *2^-4*3f/4- s i*4|0) , (A6) 

where we have set the ground state energy to eo = 0. 

The expectation values on the right-hand sides of 
Eqs. (|A5|) and (|A6|) are best evaluated in the eigenbasis 
of the Hamiltonian H . Let |a) be the set of eigenstates 
of H with energy eigenvalues e a , 



H\a) = e a \a) 



(A7) 



The eigenstates |a) = \n,m) are classified by the princi- 
pal (n) and angular momentum (m) quantum numbers 
such that in polar coordinates, the eigenfunctions 



1 



lr,Va = 



obey the following eigenvalue equation: 



1 d 



d 



m 
2^ 



rjr 



2rdr\drJ ' 2r 2 ' 2 

In terms of the matrix elements 

A a p = {a\ip\P) = (a\x + iy 



the correlation functions are 



S2l(e Q -eo) 



\A 



of) 



(A8) 

(A9) 
(A10) 

(All) 



and 



G 2 (si,s 2 ,s 3 , S4) = J2 e- K Al Q A a (3A*pA^ , (A12) 

for the case si > s 2 > S3 > S4, where we have defined 

K = e (s4 - «i) + e 7 (s3 - S4) 

+ £/j(s 2 - S3) + e Q (si - s 2 ). (A13) 



For different ordcrings of s 1 , s 2 , S3 , and S4 , similar expres- 
sions can be obtained. Although there are, in general, 
4! = 24 cases to consider, by noticing that the expecta- 
tion value for G 2 (si, s 2 , S3, S4) remains invariant under 
the exchange of s\ S3, or s 2 <^ S4, and also under the 
simultaneous exchange of si«=^s 2 and S3<=^S4, we realize 
that it is sufficient to compute G 2 (si, s 2 , S3, S4) in just 
three cases: si > s 2 > S3 > S4, si > s 2 > S4 > S3, and 
si > S3 > s 2 > S4. The remaining 21 expressions can be 
obtained from these using symmetry considerations. 

Solving the Schroedinger equation (|A9|) numerically 
and evaluating the matrix elements given by Eq. (|A10|) 
yields the correlation functions in Eqs. (|A5|) and (|A6[) . It 
should be noted that the sums in Eqs. (jA"TTT) and (|A1~2|) 
contain only a finite number of terms because of the selec- 
tion rule, A a p oc 5 mojTn/3+ i, and the fact that the bra-ket 
states with a very large difference in the respective values 
of n give negligible matrix elements due to very different 
nodal structure. 



Appendix B: Quasicondensate limit 

In this appendix, we show that the classical field ap- 
proximation correctly predicts the correlation functions 
in Eqs. and (|2"7|). in the limit -q > 1. In this limit, the 
wavefunctions of the lowest-lying states differ from zero 
significantly only for r ~ tq = yfrj, so that the Hamil- 
tonian becomes separable into the azimuthal and radial 
degrees of freedom. This has two consequences on the 
classical field calculations. 

First, the wave functions 4*™{ r ) are approximately in- 
dependent of m, while the azimuthal kinetic energy is 
reduced to 



H 



1 d 2 m 2 
'^dd 2 = 2T 2 ' 



(Bl) 



Accordingly, for the matrix elements (|A10[) we obtain 



Aaa' = S mi m'+l{ ( t>n\ r \ ( t > n 



r 5. 



m,m'+l\9n \Yn> 



(B2) 



This allows one to restrict summations in Eqs. (jAllj) and 
(|A12[) to just the leading term with n a = rip = n~ ( = 0. 

Second, the fact that the energy eigenvalues sepa- 
rate into m-independent and angular parts, 



9r 2 ' 
z 'o 



(B3) 



allows one to calculate the exponentially decaying terms 
for Gi(si,s 2 ) in Eq. (|A11[) . More precisely, the only 
relevant energy differences are 



2r 2 



1 



and 



f 2 - f 1 - — - — - 3C? 1 - f") 
~ 2-q ~ 



2r 2 



(B4) 



(B5) 
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Using Eqs. (|B2 j) -(|B5 j) . we then find that Eq. (fXTT|) re- 
duces to Gi(si, s 2 ) = ?ye^l Sl_S2 l/ 21 '. Going back to natu- 
ral units, using Eqs. (|40]) and (|42|) . this gives the quasi- 
condensate equation of state p ~ /<7, and we recover 
Eq. (|25|) of the main text. 

Considering now Eq. (|A12[) . together with the other 
required cases for time ordering, a similar albeit more 
lengthy calculation shows that Eq. (|A12[) reduces to 
Eq. ([57]) of the main text. 



Appendix C: Ideal Bose gas limit 

The limit r) <C —1 corresponds to the highly degener- 
ate ideal Bose-gas regime. In this case, the classical field 
problem can be mapped onto a two-dimensional quan- 
tum harmonic oscillator. Here we show how the classical 
field approximation recovers the correlation functions ex- 
pected for the ideal Bose gas. 

For 77 — 1, the Hamiltonian (|46[) becomes quadratic, 

H^l(p 2 x+P l) + \v\(x 2 + y 2 ), (Ci) 

and its matrix elements can be obtained from the stan- 
dard results for the quantum harmonic oscillator with 
frequency uj = ^2\rj\. We thus have 

(a\x\0) = (a\y\0) = (2a,)" 1 / 2 , (C2) 



where n a = and m a = 1 corresponds to the first excited 
state with energy e a — eo = oj = y/^rj- Then, Eq. (jAllj) 
becomes 

Gi(si, 82) = e-^- s ^-^\(a\x + iy\0)\ 2 

VWi\ 

Going back to natural units, using Eqs. (|4T)]) and 
we have p = \f^\r\\ = ^/mk 2 B T 2 /2h 2 \p\ and therefore 

Gl(z-L,Z 2 ) ^ pe -\^-^\mk B T/h 2 p = pe -\ Zl -z 2 \/l^ ^ 

which is the result for a highly degenerate ideal Bose gas. 

The calculation of the G 2 function is more elaborate 
as there are different terms on the right-hand side of 
Eq. (|A12[) to compute, for different orderings of si, s 2 , S3, 
and S4. However, the analogy with a simple harmonic os- 
cillator makes the calculation possible, leading to the re- 
covery of the Wick's theorem (valid for quadratic Hamil- 
tonians) and therefore 

G 2 (s 1 ,s 2 ,s 3 ,s 4 ) = Gi(si - s 2 )G 1 (s 3 - s 4 ) 

+ Gi(« 1 -8 4 )G 1 (« 2 -« 3 ). (C5) 

This immediately leads to the first two terms on the right- 
hand side of Eq. (fTT|) . The last (regular) term in Eq. 
([lip is identically zero in a noninteracting gas, whereas 
the third (delta-function) term, which comes from the 
commutator [ifj(s 2 ), ip*{s3)] = S(s 2 — S3), has a negligi- 
ble contribution in the highly degenerate ideal Bose-gas 
regime considered here. 
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